options(scipen = 100)

library(tidyverse)
library(data.table)
library(sf)


od <- fread("OD_MODE.TXT")
od_am <- fread("AM_OD_MODE.TXT")
od_pm <- fread("PM_OD_MODE.TXT")
colnames(od) <- c("origin_zn", "origin_code", "origin_node", 
                  "des_zn", "des_code", "des_node",
                  "Pedestrian", "Passenger_car", "Carpooled", 
                  "City_bus", "Regional_bus", "Town_bus", "Intercity_bus", "Express_bus", "Other_buses",
                  "Subway", "Rail", "Speed_rail",
                  "Taxi", "Small_freight", "Large_freight",
                  "Motorbike", "Bicycle",
                  "Other modes")

colnames(od_am) <- c("origin_zn", "origin_code", "origin_node", 
                     "des_zn", "des_code", "des_node",
                     "pedestrian", "passenger_car", "carpooled", 
                     "City_bus", "Regional_bus", "Town_bus", "Intercity_bus", "Express_bus", "Other_buses",
                     "Subway", "Rail", "Speed_rail",
                     "Taxi", "Small_freight", "Large_freight",
                     "Motorbike", "Bicycle",
                     "Other modes")

colnames(od_pm) <- c("origin_zn", "origin_code", "origin_node", 
                     "des_zn", "des_code", "des_node",
                     "pedestrian", "passenger_car", "carpooled", 
                     "City_bus", "Regional_bus", "Town_bus", "Intercity_bus", "Express_bus", "Other_buses",
                     "Subway", "Rail", "Speed_rail",
                     "Taxi", "Small_freight", "Large_freight",
                     "Motorbike", "Bicycle",
                     "Other modes")

## Code update: Hyehwa dong 1101065 -> 1101073
od %>% 
  mutate(origin_code = case_when(origin_code == 1101065 ~ 1101073,
                                 TRUE ~ as.numeric(origin_code)),
         des_code = case_when(des_code == 1101065 ~ 1101073,
                              TRUE ~ as.numeric(des_code))) %>% 
  round(digits = 0) -> od

od_am %>% 
  mutate(origin_code = case_when(origin_code == 1101065 ~ 1101073,
                                 TRUE ~ as.numeric(origin_code)),
         des_code = case_when(des_code == 1101065 ~ 1101073,
                              TRUE ~ as.numeric(des_code))) %>% 
  round(digits = 0) ->od_am

od_pm %>% 
  mutate(origin_code = case_when(origin_code == 1101065 ~ 1101073,
                                 TRUE ~ as.numeric(origin_code)),
         des_code = case_when(des_code == 1101065 ~ 1101073,
                              TRUE ~ as.numeric(des_code))) %>% 
  round(digits = 0) -> od_pm


###################################
#--admin code for our study area--#
###################################
# Jongno: 1101072, 1101054, 1101060, 1101073, 1101053, 1101061, 1101063, 1101064
#   Jung: 1102052, 1102055, 1102060, 1102054, 1102057, 1102059, 1102058

codes <- c(1101072, 1101054, 1101060, 1101073, 1101053, 1101061, 1101063, 1101064,
           1102052, 1102055, 1102060, 1102054, 1102057, 1102059, 1102058)

####################
# 1. AADT: all day
##################
#In
## I have to make sure that inflows are divided into 3 categories:
## Clean up
### 1) Outside -> CBD 
od %>% 
  select(des_code, 7:24) %>% 
  filter(des_code %in% codes) %>% 
  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_in_tot

### 2) CBD -> CBD 
od %>% 
  filter(origin_code %in% codes, des_code %in% codes) %>% 
  select(des_code, 7:24) %>% # self
  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_in_CBD

### 3) Myself -> Myself
od %>%
  filter(origin_code == des_code) %>% 
  filter(des_code %in% codes) %>% 
  select(des_code, 7:24) %>% # self
  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_in_self

## Finalise
od_in_tot[,2:19] <- od_in_tot[,2:19] - od_in_CBD[,2:19]  #1) Outside -> CBD 
od_in_CBD[,2:19] <- od_in_CBD[,2:19] - od_in_self[,2:19] #2) CBD -> CBD 

od_in_tot[od_in_tot == 0] <- NA
od_in_CBD[od_in_CBD == 0] <- NA
od_in_self[od_in_self == 0] <- NA



######################
##- With Shapefile -##
######################
cbd_rd <- read_sf("CBD_road.shp")
cbd <- st_read("CBD_boundary.shp") %>% 
  select(adm_dr_cd, adm_dr_nm_) %>% 
  rename(Code = adm_dr_cd,
         Name = adm_dr_nm_)

# Inflow: OD from Outside towards CBD
od_in_tot %>% 
  rename(Code = des_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_in_tot_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_in_tot_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_in_tot

# Inflow: OD within CBD  
od_in_CBD %>% 
  rename(Code = des_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_in_CBD_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_in_CBD_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_in_CBD

# Inflow: OD within self  
od_in_self %>% 
  rename(Code = des_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_in_self_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_in_self_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_in_self





library(classInt)
# get quantile breaks. Add .00001 offset to catch the lowest value

## Inflow: Total
breaks_qt <- classIntervals(c(min(cbd_in_tot$Count), cbd_in_tot$Count), n = 5, style = "jenks")
cbd_in_tot <- mutate(cbd_in_tot, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: CBD
breaks_qt <- classIntervals(c(min(cbd_in_CBD$Count), cbd_in_CBD$Count), n = 5, style = "jenks")
cbd_in_CBD <- mutate(cbd_in_CBD, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: 
breaks_qt <- classIntervals(c(min(cbd_in_self$Count), cbd_in_self$Count), n = 5, style = "jenks")
cbd_in_self <- mutate(cbd_in_self, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))




cbd_in_tot %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_in_tot


cbd_in_CBD %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_in_CBD


cbd_in_self %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_in_self


#ggsave("Commute_In_Outside.png", plot_in_tot, width = 12, height = 7, dpi = 200)
#ggsave("Commute_In_CBD.png", plot_in_CBD, width = 12, height = 7, dpi = 200)
#ggsave("Commute_In_Self.png", plot_in_self, width = 12, height = 7, dpi = 200)



#################################################################################
## I have to make sure that inflows are divided into 3 categories:
## Clean up
### 1) Outside -> CBD 
od %>% 
  select(origin_code, 7:24) %>% 
  filter(origin_code %in% codes) %>% 
  group_by(origin_code) %>% 
  summarise_all(funs(sum)) -> od_out_tot

### 2) CBD -> CBD 
od %>% 
  filter(origin_code %in% codes, des_code %in% codes) %>% 
  select(origin_code, 7:24) %>% # self
  group_by(origin_code) %>% 
  summarise_all(funs(sum)) -> od_out_CBD

### 3) Myself -> Myself
od %>%
  filter(origin_code == des_code) %>% 
  filter(origin_code %in% codes) %>% 
  select(origin_code, 7:24) %>% # self
  group_by(origin_code) %>% 
  summarise_all(funs(sum)) -> od_out_self

## Finalise
od_out_tot[,2:19] <- od_out_tot[,2:19] - od_out_CBD[,2:19]  #1) Outside -> CBD 
od_out_CBD[,2:19] <- od_out_CBD[,2:19] - od_out_self[,2:19] #2) CBD -> CBD 

od_out_tot[od_out_tot == 0] <- NA
od_out_CBD[od_out_CBD == 0] <- NA
od_out_self[od_out_self == 0] <- NA


######################
##- With Shapefile -##
######################
# Outflow: OD from Outside towards CBD
od_out_tot %>% 
  rename(Code = origin_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_out_tot_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_out_tot_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_out_tot

# Inflow: OD within CBD  
od_out_CBD %>% 
  rename(Code = origin_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_out_CBD_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_out_CBD_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_out_CBD

# Outflow: OD within self  
od_out_self %>% 
  rename(Code = origin_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_out_self_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_out_self_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_out_self

library(classInt)
# get quantile breaks. Add .00001 offset to catch the lowest value

## Inflow: Total
breaks_qt <- classIntervals(c(min(cbd_out_tot$Count), cbd_out_tot$Count), n = 5, style = "jenks")
cbd_out_tot <- mutate(cbd_out_tot, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: CBD
breaks_qt <- classIntervals(c(min(cbd_out_CBD$Count), cbd_out_CBD$Count), n = 5, style = "jenks")
cbd_out_CBD <- mutate(cbd_out_CBD, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: 
breaks_qt <- classIntervals(c(min(cbd_out_self$Count), cbd_out_self$Count), n = 5, style = "jenks")
cbd_out_self <- mutate(cbd_out_self, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))



cbd_out_tot %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_out_tot


cbd_out_CBD %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_out_CBD


cbd_out_self %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_out_self


#ggsave("Commute_Out_Outside.png", plot_out_tot, width = 12, height = 7, dpi = 200)
#ggsave("Commute_Out_CBD.png", plot_out_CBD, width = 12, height = 7, dpi = 200)
#ggsave("Commute_Out_Self.png", plot_out_self, width = 12, height = 7, dpi = 200)


#####################################################################
## I have to make sure that inflows are divided into 3 categories:
## Clean up
### 1) Outside -> CBD 
od_am %>% 
  select(des_code, 7:24) %>% 
  filter(des_code %in% codes) %>% 
  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_am_in_tot

### 2) CBD -> CBD 
od_am %>% 
  filter(origin_code %in% codes, des_code %in% codes) %>% 
  select(des_code, 7:24) %>% # self
  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_am_in_CBD

### 3) Myself -> Myself
od_am %>%
  filter(origin_code == des_code) %>% 
  filter(des_code %in% codes) %>% 
  select(des_code, 7:24) %>% # self
  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_am_in_self

## Finalise
od_am_in_tot[,2:19] <- od_am_in_tot[,2:19] - od_am_in_CBD[,2:19]  #1) Outside -> CBD 
od_am_in_CBD[,2:19] <- od_am_in_CBD[,2:19] - od_am_in_self[,2:19] #2) CBD -> CBD 

od_am_in_tot[od_am_in_tot == 0] <- NA
od_am_in_CBD[od_am_in_CBD == 0] <- NA
od_am_in_self[od_am_in_self == 0] <- NA


######################
##- With Shapefile -##
######################
# Outflow: OD from Outside towards CBD
od_am_in_tot %>% 
  rename(Code = des_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_am_in_tot_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_am_in_tot_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_am_in_tot

# Inflow: OD within CBD  
od_am_in_CBD %>% 
  rename(Code = des_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_am_in_CBD_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_am_in_CBD_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_am_in_CBD

# Outflow: OD within self  
od_am_in_self %>% 
  rename(Code = des_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_am_in_self_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_am_in_self_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_am_in_self

library(classInt)
# get quantile breaks. Add .00001 offset to catch the lowest value

## Inflow: Total
breaks_qt <- classIntervals(c(min(cbd_am_in_tot$Count), cbd_am_in_tot$Count), n = 5, style = "jenks")
cbd_am_in_tot <- mutate(cbd_am_in_tot, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: CBD
breaks_qt <- classIntervals(c(min(cbd_am_in_CBD$Count), cbd_am_in_CBD$Count), n = 5, style = "jenks")
cbd_am_in_CBD <- mutate(cbd_am_in_CBD, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: 
breaks_qt <- classIntervals(c(min(cbd_am_in_self$Count), cbd_am_in_self$Count), n = 5, style = "jenks")
cbd_am_in_self <- mutate(cbd_am_in_self, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))



cbd_am_in_tot %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_am_in_tot


cbd_am_in_CBD %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_am_in_CBD


cbd_am_in_self %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_am_in_self


#ggsave("Commute_AM_Outside.png", plot_am_in_tot, width = 12, height = 7, dpi = 200)
#ggsave("Commute_AM_CBD.png", plot_am_in_CBD, width = 12, height = 7, dpi = 200)
#ggsave("Commute_AM_Self.png", plot_am_in_self, width = 12, height = 7, dpi = 200)

############################################################################
## I have to make sure that inflows are divided into 3 categories:
## Clean up
### 1) Outside -> CBD 
od_pm %>% 
  select(origin_code, 7:24) %>% 
  filter(origin_code %in% codes) %>% 
  group_by(origin_code) %>% 
  summarise_all(funs(sum)) -> od_pm_out_tot

### 2) CBD -> CBD 
od_pm %>% 
  filter(origin_code %in% codes, des_code %in% codes) %>% 
  select(origin_code, 7:24) %>% # self
  group_by(origin_code) %>% 
  summarise_all(funs(sum)) -> od_pm_out_CBD

### 3) Myself -> Myself
od_pm %>%
  filter(origin_code == des_code) %>% 
  filter(origin_code %in% codes) %>% 
  select(origin_code, 7:24) %>% # self
  group_by(origin_code) %>% 
  summarise_all(funs(sum)) -> od_pm_out_self

## Finalise
od_pm_out_tot[,2:19] <- od_pm_out_tot[,2:19] - od_pm_out_CBD[,2:19]  #1) Outside -> CBD 
od_pm_out_CBD[,2:19] <- od_pm_out_CBD[,2:19] - od_pm_out_self[,2:19] #2) CBD -> CBD 

od_pm_out_tot[od_pm_out_tot == 0] <- NA
od_pm_out_CBD[od_pm_out_CBD == 0] <- NA
od_pm_out_self[od_pm_out_self == 0] <- NA


######################
##- With Shapefile -##
######################
# Outflow: OD from Outside towards CBD
od_pm_out_tot %>% 
  rename(Code = origin_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_pm_out_tot_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_pm_out_tot_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_pm_out_tot

# Inflow: OD within CBD  
od_pm_out_CBD %>% 
  rename(Code = origin_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_pm_out_CBD_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_pm_out_CBD_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_pm_out_CBD

# Outflow: OD within self  
od_pm_out_self %>% 
  rename(Code = origin_code) %>% 
  reshape2::melt(id = c("Code"), 
                 variable.name = "Mode", value.name = "Count") -> cbd_pm_out_self_long

cbd %>% 
  mutate(Code = as.integer(as.character(Code))) %>% 
  left_join(cbd_pm_out_self_long, by = "Code") %>%
  filter(Code %in% codes) %>% 
  na.omit() -> cbd_pm_out_self

library(classInt)
# get quantile breaks. Add .00001 offset to catch the lowest value

## Inflow: Total
breaks_qt <- classIntervals(c(min(cbd_pm_out_tot$Count), cbd_pm_out_tot$Count), n = 5, style = "jenks")
cbd_pm_out_tot <- mutate(cbd_pm_out_tot, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: CBD
breaks_qt <- classIntervals(c(min(cbd_pm_out_CBD$Count), cbd_pm_out_CBD$Count), n = 5, style = "jenks")
cbd_pm_out_CBD <- mutate(cbd_pm_out_CBD, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))
## Inflow: 
breaks_qt <- classIntervals(c(min(cbd_pm_out_self$Count), cbd_pm_out_self$Count), n = 5, style = "jenks")
cbd_pm_out_self <- mutate(cbd_pm_out_self, Breaks = cut(Count, breaks_qt$brks, include.lowest=T,dig.lab=8))



cbd_pm_out_tot %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_pm_out_tot


cbd_pm_out_CBD %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_pm_out_CBD


cbd_pm_out_self %>% 
  ggplot() +
  geom_sf(aes(fill = Breaks)) +
  #geom_sf_text(aes(label = Name), colour = "black", size = 1) +
  facet_wrap(~Mode) +
  scale_fill_brewer(palette = "OrRd")  +
  theme_bw() +
  theme(#legend.position = "bottom",
    legend.title=element_text(size=12),
    legend.text=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.text.x=element_text(size=10, angle=90),
    axis.title=element_text(size=15),
    strip.text = element_text(size = 11)
  ) -> plot_pm_out_self


#ggsave("Commute_PM_Outside.png", plot_pm_out_tot, width = 12, height = 7, dpi = 200)
#ggsave("Commute_PM_CBD.png", plot_pm_out_CBD, width = 12, height = 7, dpi = 200)
#ggsave("Commute_PM_Self.png", plot_pm_out_self, width = 12, height = 7, dpi = 200)
































od %>% 
  select(origin_code, 7:24) %>% 
  filter(origin_code %in% codes) %>% 
  #  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_out_tot

### 2) CBD -> CBD 
od %>% 
  filter(origin_code %in% codes, des_code %in% codes) %>% 
  select(origin_code, 7:24) %>% # self
  #  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_out_CBD

### 3) Myself -> Myself
od %>%
  filter(origin_code == des_code) %>% 
  filter(origin_code %in% codes) %>% 
  select(origin_code, 7:24) %>% # self
  #  group_by(des_code) %>% 
  summarise_all(funs(sum)) -> od_out_self




table(od_out_CBD[,-1])
round(prop.table(od_out_CBD[,-1]) * 100 ,2)



od_out_CBD
od_out_self

